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ABSTRACT 

We search for evidence of dark matter in the Milky Way by utilizing the stellar number density 
distribution and kinematics measured by the Sloan Digital Sky Survey (SDSS) to heliocentric distances 
exceeding ^10 kpc. We employ the cylindrically symmetric form of Jeans equations and focus on the 
morphology of the resulting acceleration maps, rather than the normalization of the total mass as 
done in previous, mostly local, studies. Jeans equations are first applied to a mock catalog based 
on a cosmologically derived iV-body + SPH simulation, and the known acceleration (gradient of 
gravitational potential) is successfully recovered. The same simulation is also used to quantify the 
impact of dark matter on the total acceleration. We use Galfast, a code designed to quantitatively 
reproduce SDSS measurements and selection effects, to generate a synthetic stellar catalog. We apply 
Jeans equations to this catalog and produce two-dimensional maps of stellar acceleration. These 
maps reveal that in a Newtonian framework, the implied gravitational potential cannot be explained 
by visible matter alone. The acceleration experienced by stars at galactocentric distances of ~20 
kpc is three times larger than what can be explained by purely visible matter. The application of 
an analytic method for estimating the dark matter halo axis ratio to SDSS data implies an oblate 
halo with qoM = 0.47 ± 0.14 within the same distance range. These techniques can be used to map 
the dark matter halo to much larger distances from the Galactic center using upcoming deep optical 
surveys, such as LSST. 

Subject headings: stars: kinematics and dynamics — stars: statistics — Galaxy: general — Galaxy: 
kinematics and dynamics — Galaxy: structure — Galaxy: halo 



1. INTRODUCTION 

Determining the dark matter content of the Milky Way 
has important implications for fields ranging from the- 
ories of galaxy formation and evolution to experimen- 
tal physics. Lately, there has been renewed interest in 
an old concept - applying Jeans equations to stars in 
the Milky Way to infer the underlying mass distribution 
(|Jeanslll915l : lOortH T932L This technique statistically es- 
timates the gravitational potential using observable stel- 
lar kinematics, rather than accelerations that are rarely 
detectable. 

Jeans equations follow from the collisionless 
Boltzmann equatio n : for a detailed derivation see 
iBinney fe Tremainel (|1987l ). Using cylindrical coordi- 
nates and assuming an axi-symmetric and steady-state 
system, the accelerations in the radial (R) and vertical 
(Z) directions can be expressed in terms of observable 
quantities: the stellar number density distribution, v, 
the mean azimuthal (rotational) velocity v^, and four 
velocity dispersions, a^, o\r.r, crzz, and <trz (ah as 
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Given accelerations clr{R,Z) and az{R, Z), i.e. the 
gradient of the gravitational potential, the dark matter 
contribution can be estimated after accounting for the 
contribution from visible matter. 

Traditionally, such studies were limited by data to 
the s o lar neighborhood (within ~ 150 pc, e.g. iKaptevnl 
[1921 lOortl fl96Ct IBahcalU [l98l . The main conclu- 
sion drawn from such local studies is that dark mat- 
ter contributes a small (of the order 10%) but signifi- 
cantly detected fraction of mass in the solar neighbor- 
hood (corre sponding to about 0.01 M pc~ 3 , or 0.38 
GeV cm" 3 ; IKuiiken fe Gilmordll98i ICreze et al.l H99i 
IHolmberg fc Flvnnll2000D ~ 

Several groups have extended the se studies to a few 
kilop ar sec from the p lane of t he disk (IKuiiken fc Gilmon 
1991 ISiebert et alJ [2001 IHolmberg fc Flvnnl 1200 
Smith et all l2012t IBovv et all 120121 ). Recently 



Garbari et al.l (|2012l ) used a sample of 2000 K dwarf 



stars that extend to 1 kpc above the plane of the disk 
and estimated the local dark matter density distribu- 
tion pdm = (0.022 ± 0.015) M Q pc~ 3 . Using kinematic 
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data for ~400 thick disk stars at distances of a f ew kp c 
from the Galactic plane from iMoni Bidin et all (|2012D , 
IBovv fe Tremaind (j2012l ) estimated p dm = (0.008 ± 
0.002) M Q pc- 3 . 

It has been difficult to extend these measurements to 
distances beyond a few kilopa rsec from the solar neigh- 
borhood (jvan der M arel 199l|). Recently, using a sample 
of ~2, 500 blue h orizontal branch stars from SDSS DR6, 
iXue et al.1 ()2008l ) found an estimate of the Milky Way's 
circular velocity curve at ^60 kpc that implied the ex- 
istance of dark matter. Using a spherical approxima- 
tion of Jea n s equ ations and extending the analysi s of the 
IXue et all (|2008l ) sample, iSamurovii fc Lalovicl (f20TTh 
also concluded that the Newtonian model without dark 
matter cannot fit the observed velocity dispersion profile. 
They also tested various MOND models and concluded 
that these fit the data as well. 

Here we extend these studies and introduce a 
novel multi-dimensional application of Jeans equations 
made possible by the Sloan Digital Sky SurvejO data 
(jYork et al.1 l2000i hereafter SDSS). Due to substantial 
SDSS sky coverage and accurate multi-color photom- 
etry to faint limits, the stellar number density distri- 
bution and stellar kinematics were mapped out using 
numerous main sequence stars detected to galactocen- 
tric di s tances of ~20 kp c (jJuric et al.H2008t llvezic et"aLl 
I2008al iBond et all 12010 ) . The extent of these maps is 
sufficiently large that it is possible to investigate stel- 
lar acceleration via Jeans equations. Most importantly, 
while the spatial derivatives of the velocity dispersion 
are extremely difficult to constrain from the local solar 
neighborhood, they can be directly measured using SDSS 
data. We discuss here the following main questions: 

• How does the inclusion or exclusion of a dark mat- 
ter component affect the morphology of the stellar 
acceleration maps? 

• Is it possible to recover the known accelerations by 
applying Jeans equations to a realistic simulated 
galaxy that has a merger history and is not per- 
fectly axi-symmetric? 

• Are stellar acceleration maps derived from SDSS 
data consistent with expectations based only on 
visible matter? 

This paper provides a brief summary of our analysis; 
a detailed discussion will be presented elsewhere (Loeb- 
man et al. , in prep). In Sj^l wc describe the simulation 
employed in this work and answer the first two questions. 
The main result of this work, an answer to the third ques- 
tion, is presented in JJJ We summarize and discuss our 
results in S]4j 

2. TESTING THE METHODOLOGY 

2.1. N-body+SPH Simulation 

To test the Jeans equation approach, we apply our 
analysis tools to a simulation with known accelera- 
tions and velocities. We use a cosmologically derived 
(WMAP3. fSpergel et aT1l2003f) Milky Way-mass galaxy 
evolved for 13.7 Gyr using the parallel 7V-body+SPH 
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Figure 1. A comparison of the acceleration in the Z direction 
when all contributions are included (star, gas, and dark matter 
particles; top panel) to the result without dark matter (middle 
panel). The acceleration is expressed in units of 2.9x 10~ 13 km/s 2 . 
The ratio of the two maps is shown in the bottom panel. The 
importance of the dark matter increases with the distance from the 
origin; at the edge of the volume probed by SDSS (R ~ 20 kpc, 
Z ~ 10 kpc), the total acceleration in the analyzed simulation is 
about 3 times larger than the contribution from the visible matter. 
The maps are limited to the volume explored by SDSS data, as 
indicated by the diagonal lines encapsulating the colored pixels. 



code GASOLINE (jWadslev et al.|[200l . which contains 
realis t ic gas, cool i ng an d stellar feedback ( Stinson et al.l 
l200l IShen et al.l l2009t IGhristensen et al.l 120121 ). We 
track the galaxy's formation and evolution using the 
volume renormaliza t ion technique dKatz fc Whitel 119931 : 
iPontzen et all 120081: iGovernato et al.l 12012ft . ~ Our simu- 
lated galaxy includes a stellar halo, which is built up 
primarily duri ng the merging process in a ACDM cos- 
molo gy (e.g. iBullock fc Johnston! r2005: Zol otov et al.1 



GASOLINE simultaneously calculates the potential 
and the acceleration that particles feel; force cal- 
culations are consistent with ot her state-of-the-art 
cosmological gas-dynam ical codes dPower et al.l 120031 : 
IScannapieco et al.H2012ft. The typica l RMS acceleration 
error is ~ 0.2% (jWadslev et al.lf2004l ). The average stel- 
lar particle mass is ~ 5800 M & and the dark matter 
particle mass is 1.3 x 10 5 M©, with a dark matter soft- 
ening length of 173 pc. At redshift of zero, the sim- 
ulated galaxy has a virial radius, defined at pj p cr it = 
100, of 227 kpc (versus th e Milky Way's 200 kpc, see 
iBovlan-Kolchin et al.ll20ilft and a virial mass of 7 x 10 11 
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M 9 (versus the Milky Way's 1.0 
IXue et al.ll200lllKlvpin et al.ll200l . of this, 7% is in gas, 
6% is in stars, and 87% is in dark matter. Dark mat- 
ter consists of 36% of the total mass in the region cor- 
responding to the solar neighborhood (7 < R/kpc < 9 
& < Z/kpc < 1). The simulated galaxy is approxi- 
mately rotationally symmetric (a total matter axis ratio 
b : a > 0.9 within 100 kpc, and a stellar matter axis ra- 
tio b : a > 0.95 at i?=10 kpc), has a i?-band disk scale 
length of ~3.1 kp c (versus the Milky Way's 3.6 kpc, see 
Uuric et ahl 12008 ) and corresponding bulge to disk ratio 
of 0.33 (| Jonssonl l2006f ) . and the circular velocity at 2.2 
disk scale len gths is ^208 km/ s (versus the Milky Way's 
220 km/s, see lXue et aLll2008l) . We draw direct compar- 
isons between the simulation and the Milky Way as these 
structural parameters are within 15% of one another. 

For the sake of brevity, we focus our presentation on 
the acceleration in the Z direction, az] detailed analysis 
of the acceleration in the R direction, a^, will be pre- 
sented in a subsequent paper. The top panel of Figure [1] 
shows az within the simulation. For comparison, the 
middle panel shows an analogous map when only contri- 
butions from star and gas particles are included. As is 
evident, there are substantial differences in the morphol- 
ogy of the two maps; the bottom panel demonstrates that 
the effect of dark matter on the resulting acceleration in- 
creases quickly towards the outer parts of the galaxy; for 
example, the ratio of accelerations is doubled by R=8 
kpc and Z=6 kpc. These distances are probed by SDSS 
- hence our results suggest that the effect of dark matter 
on stellar acceleration may be uncovered in SDSS data. 

2.2. Application of Jeans Equations to the Simulation 

Here we verify that the known gravitational acceler- 
ations in the simulation can be recovered using Jeans 
equations. We bin the star particles into 1 kpc by 1 kpc 
rectilinear pixels in R — Z space and limit our analy- 
sis to bins containing at least 100 star particles. Using 
weights proportional to the mass of each star particle, 
we create maps of ln(z/) and kinematic quantities. The 
spatial derivatives of these quantities at the position of 
each pixel are computed by first fitting a second order 
polynomial in R and Z using the eight neighboring pix- 
els, and then taking derivatives analytically; edge pixels 
are discarded to minimize the impact of fitting errors. 

We compare the acceleration in the Z direction com- 
puted using eqs. [T] and [5] to the true acceleration in Fig- 
ure [2j The recovered acceleration map is close to the 
true map: the distribution of pixel values from the bot- 
tom panel is centered on 1.05 with a root-mean-square 
scatter (rms) of 0.24. Similar agreement is obtained for 
the acceleration in R direction, which is centered on 1.02 
with a rms value of 0.38. 

We note that Jeans equations recover meaningful ac- 
celeration maps even though the simulation is neither 
perfectly rotationally symmetric nor steady state. Given 
this precedence, we apply the same technique to syn- 
thetic SDSS data. 

3. RESULTS 

A direct application of Jeans equations to SDSS data 
would be complicated because of the complex selection 
effects and the impact of substructure. Instead, we em- 
ploy a catalog generated using the code Galfast, which is 




Figure 2. Top panel: the acceleration in the Z direction derived 
using cq. [2]and expressed in units of 2.9 X 10 -13 km/s 2 . Bottom 
panel: the ratio of the top panel to the true acceleration in the Z 
direction (top panel of Figure [TJ. 
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Figure 3. The three observables as implied by SDSS data (num- 
ber density distribution in the top panel, and the two velocity 
dispersions in the other panels, as marked, both expressed in units 
of 2 X 10 4 (km/s) 2 ) that are required to compute az(R, Z) using 
eq. [2] The pixels at Z < 1 kpc are unreliable due to SDSS satu- 
ration at r = 14; additionally we restrict our fits to regions with 8 
adjacent pixels. 
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Figure 4. The results of applying Jeans equations to the SDSS 
observations simulated using Galfast. The top panel shows a map 
of the acceleration in the Z direction expressed in units of 2.9 X 
10~ 13 km/s 2 . The middle and bottom panels show the ratio of the 
map from the top panel and the two model-based maps shown in 
the top two panels in Figure [T] 

designed to quantitatively reproduce the SDSS measure- 
ments, volume coverage, selection and other instrumen- 
tal effects. It incorporates smooth models for the stellar 
number density distribution, metallicity distribution and 
kinematics based on tomographic analysis of SPSS dat a 
(jJuric et all 120081: llvezic dTl51l2008at iBond et alJ[20loh . 
The most important ingredients for this st udy are the 
best-fi t power-law halo given by eq. 24 in lJuric et al.l 
((2001 . with power-law index n H = 2.77 ± 0.2, axis ratio 
qh = 0.64 ±0.1, and halo velo city ellipsoid that is in- 
variant in spherical coordinates ([Bond et alj |2010). with 
a r = 141 km/s, ag = 75 km/s, and a<f, = 85 km/s. While 
the same analytic models could be used to directly gen- 
erate acceleration maps, we use mock catalogs instead 
to account for shot noise due to a finite number of stars 
and edge effects and to ensure that the same analysis 
code verified on the simulation is used when processing 
real data. We emphasize that the cylindrical symmetry 
built into Galfast is fully consistent with SDSS data ( af- 
ter local substructure is masked, see lJuric et al.ll2~008D . 

Using Galfast, we generate a flux- limited catalog with 
14 < r < 21 and mimic the SDSS sky footprint by only 
considering high Galactic latitudes (\b\ > 30°). A halo- 
like sample of 0.61 million stars is selected using M r > 4 
and 0.25 < g — r < 0.35. This sample is dominated 
by low-metallicity main sequence F stars. The catalog 
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Figure 5. A pixel-by-pixcl comparison of the acceleration values 
implied by the SDSS data and the two model predictions that in- 
clude (black lines) and do not include (red lines) contributions from 
dark matter. The top panel corresponds to the az acceleration 
maps shown in Figure [4] and the bottom panel to the an accel- 
eration maps. The model-based acceleration maps that include a 
dark matter contribution provide a significantly better description 
of the acceleration maps derived from the SDSS data. 

lists true positions, absolute magnitudes, velocities and 
metallicity, as well as corresponding simulated SDSS ob- 
servations convolved with measurement errors. 

Figure [3] shows the three observables constructed us- 
ing the Galfast catalog that arc required to compute az 
using eq. [5J As expected, both velocity dispersions are 
smoothly varying, and the cr 2 RZ is at most 40% of the 
maximal a 2 zz value. 

The top panel of Figure 2] shows the resulting az map. 
The middle and bottom panels show the ratio of this 
map compared to the simulation data in the two maps 
shown in the top and middle panels in Figure [T] It is 
evident that the acceleration map implied by SDSS data 
is much closer to the model-based acceleration map that 
includes contributions from both baryons and dark mat- 
ter. When dark matter is not included, the model-based 
map has a different shape than the map derived from 
SDSS data: it under-predicts the observed acceleration 
by roughly a factor of three at R ~ 17 kpc and Z ~ 9 
kpc. A pixel-by-pixel comparison of the two ratio maps is 
shown in the the top panel in Figured] the fact that the 
shape of the observed acceleration map is better matched 
when dark matter is included is seen as a much narrower 
histogram. These results are not unique to az', the bot- 
tom panel in Figure [5] shows the same result for an, with 
even larger deviations between the observed map and 
the model-based map that does not include dark matter 
contribution. 

Our preliminary analysis of other simulations indicate 
that the observed acceleration map cannot be reproduced 
by simply increasing the amount of baryons and not in- 
cluding dark matter: the key difficulty is to reproduce 
the strong acceleration at 10 — 20 kpc from the Galactic 
center while simultaneously maintaining the shape of the 
acceleration map throughout the probed volume. 

4. DISCUSSION AND CONCLUSIONS 



The Shape of the Milky Way Dark Matter Halo 



5 



We have shown that the kinematics of stars can be 
used to probe the dark matter distribution in the Milky 
Way. To do so, we used a cosmologically derived N- 
body+SPH Milky Way-like simulation to demonstrate 
that Jeans equations are capable of recovering the un- 
derlying gravitational potential despite small deviations 
from an implied steady-state axi-symmetric system. The 
same simulation was also used to estimate the dark mat- 
ter contribution to the resulting acceleration maps. 

We also generated a synthetic SDSS catalog using Gal- 
fast, a code designed to quantitatively reproduce mea- 
surements of spatial density and kinematics of Milky Way 
stars, as well as SDSS volume coverage, selection and 
other instrumental effects. With this catalog, we created 
acceleration maps using Jeans equations. The morphol- 
ogy of these maps provides strong evidence for the pres- 
ence of dark matter. This evidence is not sensitive to 
the overall dark matter to baryon mass ratio (that is, 
the normalization of the acceleration maps), which can 
be considered as a fine-tuning model parameter, but is 
robustly derived from the shape of observed acceleration 
maps. 

While the method presented here does not yet pro- 
vide errors for the total estimate of dark matter in the 
Milky Way, it docs provide an expanded view of the dis- 
tribution of dark matter beyond the solar neighborhood. 
Significantly, it allows us to consider how the distribu- 
tion of dark matter evolves both radially and vertically in 
the context of a large, obscrvationally motivated dataset. 
Given a large number of simulations with different dark 
matter halo properties, it will be possible to constrain 
the overall dark matter to baryonic mass ratio and the 
shape of the Milky Way's dark matter halo. 

Before such a suite of simulations is generated, we can 
approximately estimate the sh ape of the dark matter halo 
using the method proposed bv Ivan der Marell(fT99l . He 
showed that for an oblate dark matter halo with a den- 
sity distribution similar to a singular isothermal sphere, 
it is possible to derive the dark matter axis ratio from 
the stellar halo axis ratio and the velocity ellipsoid (via 
application of Jeans equations; see his eqs. 1 and 21 
and Figure 1). Using SDSS measurements for the stellar 
halo axis ratio and the velocity ellipsoid (see 33, we es- 
timated that the minor to major axis ratio of the Milky 
Way's dark matter halo is qnu = 0.47±0.14. Compared 
to the stellar halo axis ratio, qn = 0.64 ± 0.1, the dark 
matter halo is more oblate. 

The iV-body+SPH simulation used in this work sup- 
ports the properties of the dark matter halo assumed in 
the van der Marel's method; in particular, its density 
decrease is consistent with expected l/[R 2 + (Z/qDAi) 2 ] 
within 30 kpc from the origin (with the power-law index 
decreasing from —2 to about —2.8 at larger distances). 
At the same time, the stellar halo is consistent with 
Uh = 2.77 measured by SDSS in the same distance range. 
While it is premature to declare qoM = 0.47 ± 0.14 as a 
robust measurement of the dark matter halo shape, it is 
encouraging that the simulation is at least qualitatively 
consistent with SDSS data in so many aspects. 

Moreover, it is possible to go beyond Jeans equations 
to use stellar kinematics t o probe the full pha se space 
distribution of stars. As IValluri et al.l (|2012l ) demon- 
strated, orbital spectral analysis can be used to deter- 
mine not only the shape of the inner halo, but the global 



shape of the Galactic halo as well. Their technique pro- 
vides a complementary tool to the method presented here 
for constraining the potential and the stellar distribution 
function. 

We note, there is a considerable range of Milky Way 
dark matter halo axis ratios presented in the literature. 
The reported axis ratio ranges from >0.7 a t the low en d 
(jlbata et al.ll200lD to 5/3 at the high end (jHelmil 120041) . 
Sometimes, triaxial mod els are incorporated, such as by 
lLaw fc Majewsk i (2010) who used observations of Sgr 
tidal stream and N-body modeling to conclude that the 
Milky Way has a triaxial dark matter halo that is nearly 
an oblate ellipsoid whose minor axis i s cont ained within 
the Galactic disk plane. iLux et all (|2012( ) used satel- 
lite galaxies to find the Milky Way dark matter halo is 
oblate, while iBanerjee fc Jogl ()2011h used observations of 
HI and claimed variation of the axis ratio, with the halo 
increasingly prolate at large distances. 

The full promise of various methods for estimating the 
geometry of dark matter halos will be reac hed by up- 
comi ng next-gener ation surveys, such a s Gaia (jPerrvmanl 
I2002T ) and LSST (jivezic et al.l l2008bO . Gaia will pro- 
vide measurements of geometric distances and kinemat- 
ics with a similar faint flux limit as SDSS, but with much 
smaller errors. LSST will obtain photometric distances 
and kinematics of comparable accuracy to those of Gaia 
at Gaia's faint limit, but extend them deeper by about 
5 mag. With LSST and Gaia, it will be possible to ex- 
tend Galactic potential studies such as the one presented 
here to distances roughly 10 times larger than those pos- 
sible with SDSS data, revolutionizing our understanding 
of the Milky Way in the process. 
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